Few-Body Systems 0, 1-20 (2008) p 

bystems 



© by Springer- Vcrlag 2008 
Printed in Austria 



Operation of Faddeev-Kernel in 
^ ; Configuration Space 

(N 

S. Ishikawa* 

Department of Physics, Science Research Center, Hosei University, 2-17-1 Fujimi, Chiyoda, 
Tokyo 102-8160, Japan 



■^j- ; Abstract. We present a practical method to solve Faddeev three-body equa- 

tions at energies above three-body breakup threshold as integral equations in 
coordinate space. This is an extension of previously used method for bound 
O ! states and scattering states below three-body breakup threshold energy. We 

show that breakup components in three-body reactions produce long-range 
effects on Faddeev integral kernels in coordinate space, and propose numeri- 
cal procedures to treat these effects. Using these techniques, we solve Faddeev 
equations for neutron-deuteron scattering to compare with benchmark solu- 
tions. 
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So far, a number of numerical methods to solve Faddeev three-body equations 
for energies above three-body breakup threshold have been developed and then 
applied to a system of nucleon-deuteron, which is considered as one of the most 
basic quantum three-body systems [TJ [2] . These methods are classified into two 
groups: either to solve coupled integral equations for scattering amplitudes in 
momentum space or to solve coupled partial differential equations for wave func- 
tions in coordinate space. 

In this paper, we will present a different approach for three-body scattering 
problem above the breakup threshold energy, in which we solve the Faddeev 
integral equations for wave functions in coordinate space. This approach has 
been successfully applied to calculations of the three-nucleon bound states [U Hj 
and low-energy three-nucleon scattering below the breakup threshold energy with 
inclusion of three-nucleon forces [4] and the long-range Coulomb interaction [5j [6]. 

Integral equations for scattering problems are generally written in the form 
of inhomogeneous linear equations. In the previous works, we applied an iterative 
method called the method of continued fraction (MCF) to solve such equations, 
whose details are given in Refs. [21 [7] and references therein. A basic procedure of 
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the algorithm in the MCF is to operate the integral kernel to a function made in 
a preceding step, as are those in most iterative methods. It is thus essential to es- 
tablish precise operations of integral kernels for solving the equations accurately, 
which is main subject of this paper. 

The existence of three-body breakup channels causes some difficulties in 
three-body calculations. In the momentum space approach, for example, the 
effects appear as logarithmic singularities and discontinuities by a step function 
in the integral kernels of the equations so that we need to perform the integration 
very carefully [8j. In the differential equation approach, due to the breakup ef- 
fects one needs to set boundary conditions at very long distance, the order of tens 
or hundreds times larger than the range of interaction potentials [9J \10\ [TTT [T2] . 
Since we treat the wave functions as solutions of the Faddeev integral equations, 
the long-range behavior should appear in the integral kernel. In the present pa- 
per, we will describe how this behavior appears in our kernel, and how to treat 
it. 

Basic notations and steps of the kernel operation in detail are explained in 
Sec. [2] for a simple three-body system. In Sec. [31 we show numerical examples of 
the kernel operation to a model function emphasizing some techniques to treat 
breakup effects, and then compare our calculations with benchmark tests [H [2] . 
Finally, we give a summary in Sec. [H 

2 Formulation 

2.1 Notations 

Let us consider a system of three identical particles (nucleons) 1, 2, and 3. We 
use sets of Jacobi coordinates {xi,yi} defined as 



where (i,j,k) denotes (1,2,3) or its cyclic permutations and rj is the position 
vector of the nucleon i. For simplicity, we assume that the nucleons j and k 
interact via a short range pair wise potential Vi = V(xi), where Xj = \xi\, and 
that the potential supports a s-wave bound state (the deuteron) of energy E^, 
whose radial part wave function is denoted by (j) d {xi). 

We are going to obtain a wave function & corresponding to a scattering pro- 
cess initiated by a state with a nucleon and a deuteron having relative momentum 
Pq. Faddeev equations for the process in the form of integral equations are 






where ^j's are Faddeev components to make \P as 



= <p l + <|> 2 + <2> 3 



(2.3) 



Si is an initial state consisting of the deuteron for the pair (j, k) and incoming 
free nucleon i, and Gi is a three-body channel Green's operator with the outgoing 
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boundary condition, 



Gi = 5 5 . (2.4) 

E + ie + **Vl + f^V 2 , -Vt 

m Xi Am Vi 1 



y a {xi,vi) = \y L (xi) ® y«(yi)]w„ > ( 2 - 7 ) 



The total energy in the three-body center of mass (cm.) frame E is given as 

3ft 2 9 3ft 2 9 , , . , 

^ = ^0 + ^ = ^-1^1, (2-5) 

where m denotes the nucleon mass. The operator P represents permutations of 
the particle numbers, 

P#i = *i+**- (2-6) 

A partial wave decomposition is performed by introducing an angular func- 
tion y Q (xi,yi), 

r\ . \ * •■ iMo 

where L denotes the relative orbital angular momentum of the pair (j, k) ; I the 
orbital angular momentum of the spectator i with respect to the cm. of the pair 
(j, k); J the total angular momentum of the three-body system (J = L + £); 
Mq the third component of Jo- The set of the quantum numbers (L,£, Jq, Mo) 
are represented by the index a. Furthermore, we use an index ao to denote an 
initial partial wave state specifically with L = 0. 

2.2 Kernel Operation 

In this subsection, we describe how to handle the operation of the Faddeev kernel 
GVP on a given function S, 

(x, y \s) = J2y«(x,y)Ux,y) (2-8) 

a 

to produce a new function <P, 

(x,y\$) = (x,y\GVP\S) 

= ^J Q (i,y)0 a (x,y), (2.9) 

a 

where we have dropped the particle number indices k) for simplicity. 

The kernel operation starts with the permutation operator P to define a 
function Xa(x, y), 

Xc{x,y) = {y a \P\S). (2.10) 
In the case of identical particles, P is nothing but a coordinate exchange operator, 



whose operations are summarized in Appendix A: 



Next step is the operation of the Green's operator G. In the case of the 
scattering problem, where E > 0, the Green's operator G possesses a pole cor- 
responding to the deuteron bound state. In order to treat this pole, we apply a 
standard subtraction method, in which we insert a trivial identity, 



l=El^«D^)(^«ol + 



i-Ei^x^oi 

CO 



(2.11) 
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between G and V in Eq. (|2.9f) . This procedure extracts an elastic contribution 
of the Green's operator [13] and leads to an expression, 



M*,V) = 5 a , ao <p d (x)^ e \y)+^> c \x,y). (2.12) 
Here, T^ e \y) represents an elastic component in the scattering given by 

T^(y)= y' 2 dy'G 0/o (y,y')J e Hy>), (2.13) 



where Goi(y, y') is a partial wave component of the free Green's operator for the 
outgoing particle, 



G Ay,y') = [y 



l 



3 £pl + ie-T e (y) 



V 



Am, 1 

-■^Poht (poy>)jt(poy<) (2.14) 



with 

3ft 2 / d 2 2 d 1(1 + 1) 
Am I dy 2 y dy y 



Ti(y) = -—1—2+-^-- • ( 2 - 15 ) 



In Eq. (|2.14p . ji(poy) is the spherical Bessel function and h^\poy) is the spher- 
ical Hankel function with the outgoing wave, where the outgoing (+) and the 
incoming (— ) spherical Hankel functions are defined with the spherical Neumann 
function n^p^y) as 

hf\x) = -n i {x)±\j l {x). (2.16) 

The function u^ e \y), which plays a role of the source for the elastic compo- 
nent in Eq. (|2.13p . is given by 

rod 

J e \y) = x 2 dx<p d (x)V(x) XQ0 (x,y). (2.17) 



The explicit expression of the Green's function Eq. (|2.14p gives the asymp- 
totic form of J-^ (y) as 

F^(y) -+ h { +\ Po y)T^, (2.18) 

y >oo 

where is the elastic T-matrix amplitude defined by 

T(e) = ~ po {S) l°°y 2d yuo(poy)^ e) (y). (2.19) 

The second term in the right hand side of Eq. (|2.12p expresses three-body 
breakup and closed-channel components in the scattering. In our formalism, these 
components are treated by expanding the Faddeev kernel with respect to a spec- 
tator particle state of momentum p, 

ui(y;p) = ^pje(py), (2.20) 
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which satisfies a complete relation 

Kv-v') = 
yy' 



dpue(y;p)ue(y';p). (2.21) 



The function <f>a' c \x,y) thereby is written as a Fourier-Bessel transformation: 

poo r , 

cp£ c) (x,y)= dpu e (y;p) Va (x;p)-5 a , ao (l ) d (x)C a ( P ) . (2.22) 
Jo 1 J 

Here, r] a (x;p) is defined as 

V a (x;p) = (x\G L \oj a ), (2.23) 
where Gl is a two-body Green's operator 

Gl = Eq + is-T L (x)-V(x) (2 ' 24) 

with 

Tl(x) = (jL + 2 A _ L ( L + 1 ) ^) . (2 . 25) 
ml a!x 2 x dx x 2 J 

The energy of the two-body subsystem £7 g is given by 

E q = E — —p 2 = —q 2 , 2.26 
4m m 

and the p-dependence of the functions arises through this relation. 

The breakup component stems from the integral of the first term in Eq. (|2.22j) 



for the range of < p < p c = <J AmE '/3h 2 '. In this range the energies of both 
the spectator particle and the two-body subsystem are positive or zero, and thus 
the integral survives at infinite values of x and y, see Eq. (I2.33P below and Refs. 
|144 115]. The rest of the integral of the first term in Eq. (|2.22p . i.e., p c < p < oo, 
as well as the second term in Eq. ()2.22j) damp for large values of x and y because 
the energy of the two-body subsystem is negative. In this sense, we call these 
components closed. 

The source term in Eq. (I2.23|) . uj a {x;p), is written as 

u a (x;p) = V(x)xa(x;p), (2.27) 

POO 

Xa(x;p)= y 2 dyui(y;p)xa(x,y). (2.28) 



The second term of the right hand side in Eq. (|2,22p appears as a counter 
part of the subtraction and C Q (p) is defined as 

C a (p) = - — / x 2 dx(t> d {x)Cj a {x-p). (2.29) 

£>q - Ed JO 

The apparent singularity in C a (p) cancels that of the two-body Green's operator 
Gl, which will be numerically shown in the following section, and thus, we can 
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apply a standard quadrature to perform the p-integration in Eq. f|2.22 j> as far as 
the both terms are treated together. 

In calculating r/ a (x;p), we transform Eq. (|2.23|) to an ordinary differential 
equation: 

[E q - T L {x) - V{x)\ r, a {x;p) = u a {x;p) (2.30) 



with boundary conditions 



h£\qx) {0<p<Pc) 



,QV '"»-oo\ h { +\i\q\x) ( Pc <p <oo) 
A treatment of the two-body Green's operator at three-body breakup re- 



gion, < p < p c will be described in Appendix B: We here only note that the 
asymptotic form of r] a (x;p) is given by 

r, a (x;p) -> h£\qx) V (&(g)|6«), (2-32) 

where q) is a two-body scattering solution with the standing wave bound- 

ary condition and K,l{q) is a scattering .fT-matrix for the two-body scattering 
(See Appendix B:[ ). 



The asymptotic form of 4>a' c \x,y) is evaluated by the saddle-point approxi- 



mation [Mj [T5] as 

€< c \x,y) -7 fl -e* ^ (^) ^7F^(0), (2.33) 

x—>oo,x/y fixed \ O / it ' 

where we introduce a hyper radius R and a hyper angle as 



/?: \x 2 + ~y 2 , (2.34) 



/3 

x = i?cos@, y = y -RsinO, (2.35) 



and Kq is given by 



K = JjpE. (2.36) 



B a (0) is the breakup amplitude defined as 



1 TTl I 



Here, the momenta q and p are given as 



/4 

g = K o cos0, p= J-K Q smG. (2.38) 

V o 
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3 Numerical Analyses and Results 

3.1 Model source function 

In this section, we present a numerical example of the kernel operation described 
in the preceding section with a model source function that is restricted to L = 
£ = Jo = state but carries a feature of the presence of three-body breakup 
channel similarly as the one used in Ref. [12) : 



e iK R 

Xaix ' y) = (R + Ro)^ (3 - 1} 

with Rq = 5 fm. 

We choose the 3 Si-component of the Malfliet-Tjon model as presented in 
Ref. [I] for the potential V{x) and the incident nucleon energy of -E_L a f> = 14.1 
MeV, which gives /\o=0.416 fm -1 , p c =0.480 fm -1 , and po=0-550 fm -1 . In nu- 
merical calculations below, mesh points for x- and y-variables in described in 



Appendix C: are used 



3.2 Elastic part 

In Fig. [H we plot the real part of the elastic source function oj( e \y), Eq. (|2.17p . 
calculated with the model function Eq. (|3.ip . As shown in this figure, uj^ e \y) 
reveals a long-range behavior, which is given by 

i \ e [pcV 
u {e) (y) oc (3.2) 

whose oscillation length is about 13 fm. 

In calculating the elastic component J-( e \y), we treat this long-range behav- 
ior by rewriting Eq. (|2.13|) as 

jr{e) {y) = _ n£o ( poy )T(y) + ij e<) (p y)T^ + j io (p y) (s(y) - s) , (3.3) 
where we have defined T(y), S(y), and S by 

In numerical integration in Eq. (|3.4p . we need to be careful for oscillational 
behaviors of the spherical Bessel and spherical Neumann functions as well as 
uj^ e \y). This is done by spline interpolation technique used in Ref. [5] taking 
into account of oscillational behavior of the integrand carefully. 

For the use of Eq. (|3.3p . one needs converged values of T(y) and S(y) for 
y — > oo. In Fig. we plot the real part of T{y) for an example. As is expected 
from the long-range behavior of oj( e \y), the convergence of T(y) becomes very 
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4x10 




y(fm) 



Figure 1. Real part of the elastic source function cj^ e '(y). 

slow. However, from the functional form of Eq. (|3,2p . we expect that the function 
T(y) behaves asymptotically as 

e Kpo-Pc)y 

T(y) - t + h (3.5) 

where to an d ^i are expansion coefficients and the coefficient to is considered as a 
converg ed value ofr( e ). The wave length evaluated by this equation is — ^— = 90 
fm, which is actually observed in Fig. [2j The fitting coefficients in Eq. (|3.5p are 
evaluated by a least square fit. To do this, the calculated values of T(y) in a 
range of 80 < y < yjff* (fm) are used. In Table [U the dependence of the result 
on yj 1 ^ is displayed. From the table, we set yj 1 ^ = 1000 fm to get a converged 
result in five digits of accuracy, which is denoted by the dashed line in Fig. [2J 
We remark that this result is obtained in spite of the fact the deviation of T(y) 
from the converged value is still about 0.5 % at y = 1000 fm, which is not shown 
in Fig. [2J 

Here, we consider a range of the variables {x, y} to be used in calculations. 
In the Faddeev equation, the function Xa(x,y) is always accompanied by the 
potential V(x), which means that we need to calculate this function within the 
range of potential, xr, for the x-variable. On the other hand, there is no restric- 
tion for the y- variable. In actual, Table Q] demonstrates that we need to calculate 
Xa(x,y) for a large value of y, i.e., 1000 fm. 

Suppose that we calculate Xa(x,y) by Eqs. (|2.8[) and (|2.10p with a func- 
tion £ a (x,y), which is obtained in a preceding iteration step, for a range of 
jo < x < xr,0 < y < yjlf*}- The formulae, Eqs. (|A.4p and (|A.5p . show that we 
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Deviation 




-0.070 III. 1 . 1 . 1 

200 400 600 

y(fm) 

Figure 2. The real part of the function T(y). The obtained converged value is shown by the 
dashed line. Dotted lines with the indices in the right hand side axis denote the deviation of 
the real part of T(y) from the converged value. 



Table 1. The real part 
of the fitting coefficient 
to. 



max /c m \ 

Vfit (fin) 


Re[to] 


532 


-6.5454 


622 


-6.5452 


983 


-6.5449 


1073 


-6.5448 


1163 


-6.5448 


1193 


-6.5448 



need to prepare the function £ a (x,y) for a range of 

0< y <\x R + \yfr (3.6) 

to perform the exchange operation to obtain Xa(x,y) for the above range. If we 
set x R =W fm and y^=1000 fm, this turns to be {0 < x < 1005 (fm),0 < y < 
507.5 (fm)}, which is rather huge. 

To facilitate numerical calculations, we limit the range of calculating Xa(x, y) 
to {0 < x < xr, < y < yj\f} by choosing the value of j/m adequately, and 
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approximate the value of Xa(%> y) for y > yM using a form of 



Xa(x,y) 



0<x<x R ,y>y M y 



(a (x) + 



2/5 



+ 



a 2 (x) 



y 



(3.7) 



where the coefficients a n (x) are determined by a least square fit to Xa(x,y) for 
y < an d for each value of x. 

With a choice o£xm = 10 fm and ?/m = 80 fm, by which the range for £ Q (x, y) 
becomes {0 < x < 85 (fm), < y < 47.5 (fm)}, we obtain the equivalent results 
for T(y) and its asymptotic value to the previously shown. This procedure 
reduces the amount of calculations considerably without loss of accuracy, and 
will be applied in the following analyses. 

Together with the function S(y) and its asymptotic value S calculated sim- 
ilarly, the elastic component !F^ e \y) is constructed using Eq. (j3.3|) . whose real 
part is plotted in Fig. [3l Note that the effect of the slow convergence in T(y) and 
S(y) functions appears as a small oscillation of the amplitude of J-^ e \y) with 



the wave length of 



2tt 
PO-Pc 



90 fm, which exists up to a large distance where the 



convergences of T{y) and S(y) are achieved. 



0.2 - 



or 



-0.2 - 




100 200 

y(fm) 



300 



Figure 3. The real part of the elastic function T^ B '(y). 



3.3 Breakup and closed channel parts 

First step in the calculation of the three-body breakup and closed channel con- 
tributions is the Fourier-Bessel transformation of Xa(x,y) with respect to the 
coordinate y, Eq. (|2.28[) . We again face the problem of slow convergence in the 
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y-integral due to the long-rangeness of Xa(x,y)- This is treated similarly with 
the calculation of the elastic component by writing Eq. (|2.28[) as 

rVM 

Xa(x;p) = / y dy'u e (y;p)x a (x,y') 
Jo 

ry 

+ lim / y' 2 dy'ut(y';p)xa(x,y'). (3.8) 

The first term is integrated numerically using the spline interpolation tech- 
nique [5]. Results for the real and the imaginary parts of the integrals with 
VM = 80 fm are shown in Fig. 2] (a) by the solid curves. The oscillational behav- 
ior of the curves indicates that the integrals do not converge yet. 




Figure 4. (a) The real and imaginary parts of Xa{x;p) for < p < 1 (fm -1 ) at x = 1 fm. 
The solid curves are the first term of Eq. (|3.8p with yu = 80 fm. The dashed curves are the 
full calculation, (b) The imaginary part and (c) the real part of Xa(x;p) around p = p c at x — 
1 fm for various values of yi n f ■ See the text for the details. 

In calculating the second term in Eq. (13.80 . we use the asymptotic form of 
Xa(x,y) given by Eq. (|3.7p . Now, we define a function l^ n \y;p) (n =0, 1, or 2) 

as 

ry pWcy' 
-.it ii \ ^ 



i {n \y;p)=J y /y'y'My'-,P) w ^, (3-9) 

and then express this for large values of y in a form of 

^(p) + 4 n) (p)^r (3.io) 
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with expansion coefficients and to be determined by a least square fit. 
The wave length of the oscillation of Eq. (|3.10j) with respect to y-variable depends 
on the momentum p as J^L . In a particular case of p = p c , where no oscillation 

V Pc 

occurs, a functional form of 

. h '{n) An) 

u yl/2+n y3/2+n v ' 



is used. The second term of Eq. f|3.8j) is thereby expressed as 

j2a n (x)b { ;\p). (3.12) 

n=0 

Since the functions j( n \y;p) depend only on yu and the total energy E, we 
may calculate them once in advance to start an iterative process in solving the 
Faddeev equations. 

The ^-coefficients in Eq. (13. lOh are obtained from a least square fit using 
values of X^ n \y;p) for a range up to y = yinf- To obtain accurate values of the 
coefficients, we need to include at least several oscillations in the range. Since 
the wave length of the oscillation becomes larger as p approaching to p c , the 
maximum value yi n f to get a converged result could become a huge number. 
This is illustrated in Figs. H] (b) and (c), where the dependence of the resultant 
Xa(x;p) on some selected values of yi n f is plotted. In the figures, we plot the 
real and imaginary parts of Xa{x;p) around p = p c = 0.480 fm -1 at x = 1 fm 
calculated by choosing y in j = 10 3 fm (dot-dashed curves), 10 4 fm (dotted curves), 
10 5 fm (dashed curves), and 4 x 10 5 fm (solid curves). One sees that even the 
value of yi n f = 10 3 fm is not enough to get a converged result. Numerically, it 
turns out that 4 x 10 5 fm may be good enough. The results with yi n f = 4 x 10 5 
fm are plotted as dotted curves in Fig. 0] (a). The oscillating behavior due to 
the small value of the integral maximum given by the solid curves disappears by 
taking into account of the long-range character of the source function Xa( x ,y)- 

Using thus obtained Xa(x',p), one calculates u a (x;p) from Eq. (|2.27p . and 
then solves the ordinary differential equation, Eq. (|2.30p . with the boundary con- 
ditions Eq. (|2.3ip to obtain rj a (x;p). The Numerov algorithm is applied for solv- 



ing this equation as in Refs. [16^15] with x-mesh points described in Appendix C: 
Fig. [5] displays the real (imaginary) part of the resultant r] a (x;p) function at 
x = 1 fm as thin solid (thin dashed) curve. The discontinuous singularities of thin 
curves at p = po = 0.550 fm" 1 correspond to the deuteron pole in the two-body 
Green's function. These singularities disappear when the term of (j) d (x)C a (p) in 
Eq. (I2.22p is subtracted, as shown as bold curves in the figure. 

In our formalism, the breakup amplitude B a {0) is obtained by two different 
ways. One way is to use Eq. (I2.37|) directory, which can be performed before 
solving the ordinary differential equation Eq. (|2.30p . After getting a solution of 
Eq. (|2.30p . the breakup amplitude is calculated from its asymptotic form Eq. 
(|2.32p as the second way. Both calculations agree each others, which assures the 
accuracy of the solutions of Eq. (|2.30p , and displayed in Fig. [6j In the inserts of 
Fig. [6j results with different values of yi n f in calculating Xa(x',p) are displayed 
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Figure 5. Thin curves are the real part (solid curve) and the imaginary part (dashed curve) of 
T] a {x;p) at x = 1 fm as functions of p. Bold curves are those after the subtraction of <f> [x)C a {jp)- 
term in Eq. fl2J22j. 



as in Fig. S] to see effects of the long-range properties of Xa(x, y) in the region of 
~ tt/2 region, where q ~ 0. 

Once the function rj a (x;p) is obtained, by performing the transformation Eq. 
(|2.22|) with the spline interpolation technique, we obtain the function (f>a' c \x, y). 
Together with the elastic component J-( e \y), we finally obtain (j)(x,y) by Eq. 

3.4 Comparison with the Benchmark solutions 

The formalism for the operation of the Faddeev kernel described in the preceding 
sections is easily extended to more realistic cases, with spin degrees of freedom, 
with three-body forces, etc. Accommodating the formalism in the MCF algo- 
rithm [7] , we are able to solve the Faddeev integral equations in coordinate 
space. To demonstrate the accuracy of our method, we performed calculations of 
the neutron-deuteron (n-d) scattering with the Malfliet-Tjon I-III potential, for 
which benchmark tests exist [HE]- The comparison are made in Tables [2] and [3] 
and Figs. [7] and 

In Tables [2] and El where we tabulate results of the s-wave phase shift pa- 
rameters for the n-d doublet and quartet states at the incident energies of 
4.0, 14.1, and 42.0 MeV, the calculations in the benchmark tests are denoted 
as Utrecht, Jiilich/NM, Bochum, LA/Iowa, and Hosei(Q). (See Ref. [1] for fur- 
ther references of these methods.) In the calculations indicated as Utrecht and 
Bochum, coupled two-dimensional integral equations in momentum space are di- 
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(deg) 

Figure 6. The breakup amplitude B a (&). The solid line shows the real part and the dashed line 
the imaginary part. The inserts show the behavior of j/i„/-dependence in calculating x a (x;p) 
function. The meaning of each curve is explained in the text. 

rectly solved by Pade approximant methods. Integral kernels of their equations 
consist of free three-body Green's operator, two-body i-matrix, and permuta- 
tions operators. The two-body £-matrix possesses a pole due to the deuteron, 
whose effect is treated by a subtraction method. The breakup effects appear as 
singularities in the three-body Green's function, see Ref. [8] for the details. Those 
indicated by Jiilich/NM and Hosei(Q) use a separable expansion for two-body 
t-matrix to reduce the dimension of integral equations to one, and then solve the 
resulting equations taking into account of singularities in the kernels by tech- 
nique of the contour deformation. In the calculations denoted as LA/Iowa, the 
Faddeev differential equations in configuration space are solved with boundary 
conditions for the elastic and the breakup regions of the wave functions. We 
noted that the boundary condition for the elastic channel used there does not 
include the small oscillation behavior found in Fig. 

A breakup amplitudes defined in Ref. [2], A(0), is related with our amplitude 
B a (G) for L = I = as 

(4\ 3/2 
-J p o K*B a (0). (3.13) 

In Figs. [7] and El the results of the breakup amplitude are compared for £x a b = 
14.1 MeV and 42.0 MeV, respectively. In the figures, our results for the real 
(imaginary) part are shown as the solid (dashed) curves, while those by Bochum 
and LA/Iowa groups [2], which are almost equivalent, are denoted by circles 
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Table 2. Comparison of the benchmark calculations [T] and the 
present calculations for neutron-deuteron spin-doublet phase shift 
parameters with the Malfliet-Tjon I-III potential. 



E Lab (MeV) 4.0 14.1 42.0 





Re(S) 


'/ 


Re(<5) 


n 


Re(<5) 


V 


Utrecht [T] 


143.7 


0.963 


106.5 


0.468 


41.9 


0.488 


Jiilich/NM [1] 


143.7 


0.952 


104.9 


0.460 


41.3 


0.501 


Bochum [T] 


143.7 


0.964 


105.5 


0.467 


41.3 


0.504 


LA/Iowa [1] 


143.7 


0.964 


105.4 


0.463 


41.2 


0.501 


Hosei(Q) U 


143.7 


0.964 


105.5 


0.465 


41.3 


0.502 


This work 


143.7 


0.964 


105.5 


0.466 


41.6 


0.498 



Table 3. Comparison of the benchmark calculations [T] and the 
present calculations for neutron-deuteron spin-quartet phase shift 
parameters with the Malfliet-Tjon I-III potential. 



E Lab (MeV) 4.0 14.1 42.0 





Re(<5) 


'/ 


Re(<5) 


V 


Re(<5) 


V 


Utrecht p] 


102.1 


1.000 


68.8 


0.978 


38.4 


0.898 


Jiilich/NM U 


101.1 


1.000 


68.5 


0.986 


37.2 


0.907 


Bochum fTJ 


101.6 


0.999 


69.0 


0.978 


37.7 


0.903 


LA/Iowa [I] 


101.5 


1.000 


68.9 


0.978 


37.8 


0.906 


Hosei(Q) P 


101.6 


1.000 


68.9 


0.978 


37.7 


0.903 


This work 


101.6 


1.000 


69.1 


0.976 


37.8 


0.889 



(triangles) . 

All of our results agree with the benchmark calculations better than 1 % 
level except about 2 % discrepancy for the rj parameter in the quartet state at 
42.0 MeV, which demonstrates the present formalism is promising in solving the 
three-body scattering problem at energies above three-body breakup threshold. 

4 Summary 

We have presented a method to operate the Faddeev integral kernel in coordi- 
nate space at energies where three-body breakup reactions take place. Effects of 
three-body breakup reactions appear as a long-range source contribution to the 
elastic component and to the breakup amplitudes at two-body sub-system hav- 
ing almost zero energy. Some numerical procedures are developed to treat these 
long-range behaviors. With a model source function and a model potential, we 
have displayed some numerical examples to verify the accuracy of our method. 

The procedure described in this paper can be used to solve the Faddeev 
equations in combination with an iterative algorithm to solve linear equations 
such as MCF. Solutions of the three- nucleon Faddeev equations are given for the 
Malfliet-Tjon I-III potential, and scattering phase shifts as well as the breakup 
amplitudes obtained from the solutions give a good agreement with the bench- 
mark solutions. Results for three- nucleon systems with realistic nucleon- nucleon 
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Faddeev-Kernel in Configuration Space 




30 60 
Hyper Angle (deg) 




30 60 
Hyper Angle (deg) 



Figure 7. The real part (solid line) and the Figure 8. The real part (solid line) and the 

imaginary part (dashed line) of the breakup imaginary part (dashed line) of the breakup 

amplitudes at _E.Lai)=14.1 MeV in the present amplitudes at i?La!>=42 MeV in the present 

calculation. Benchmark calculations in Ref. [2] calculation. Benchmark calculations in Ref. [2] 

are plotted as circle and triangle points. are plotted as circle and triangle points. 

interactions and three-nucleon interactions will be presented elsewhere. 

Since the integral kernel and hence the wave function in our formalism can 
be written as the sum of the elastic, three-body breakup, and closed channels, 
effects of each reaction mechanism can be easily drawn. Our formalism thus can 
be extended to treat a three-body model of nuclear reactions including three- 
body breakup reactions in such a way that the theory resembles conventional 
theories of reactions. 

Acknowledgement. This research was supported by the Japan Society for the Promotion of 
Science, under Grants-in-Aid for Scientific Research No. 13640300. The numerical calculations 
were supported, in part, by Research Center for Computing and Multimedia Studies, Hosei 
University, under Project No. Iab0003. 

Appendix A: Particle Exchange Operator P 

In this appendix, we summarize formulae to accomplish the particle exchange operator P in 
Eq. I|2.10p . Xa{x,y) in Eq. (|2.10|) is given as 

Xafay) = ^2xa, a '(m,y), (A.i) 

at' 

where 

xa, a '(x,y) = (y a \P\y a >ta>) 



S. Ishikawa 17 

l' e 

= J2J2^ 5a + b ' L ' Sc+d > e ' xa+Cyb+d J2 K "'( X > y ^ R (LlL>e')L ( A - 2 ) 

a — c— b,d 7i^o 



with 



^Wo - (-i)^ / ^i^M(^ 1 ) Va ( !tf ' + 1 ) 1/a 



x{ c d f 1 (LeOO| 7 0> <</00|tO) \ * 7 L , \ (A.:!) 



e f L 



I L f 



and 



\ \u ^^ tl P,{u). (A.4) 



Here, n denotes v2n+T; P 1 {u) is the Legendre polynomial; a; and ?/ are 



(A.5) 



Appendix B: Green's operator 

In this appendix, we first review two-body Green's operators and describe how to calculate Eq. 
JH3]). 

We define Green's operators for the outgoing (+) and the incoming (— ) boundary conditions 
with and without a potential as 

r (±) = 1 c R -n 

E q ±ie-T L (x)-V(x) y ' ' 

G ^ = E q ±ie-T L {xY (B - 2) 

These satisfy resolvent relations 

r.(±) _ r-(±) , ^(±)i /r .(±) _ r-(±) , ^(±) vr .(±) m 

b L 0,L ' L 0,L 0,L 0,L L • l 13 -^ 

Two-body scattering wave functions corresponding to the outgoing and the incoming bound- 
ary conditions ') satisfy the Lippmann-Schwinger equations 

\^ ) ) = \U)+G^V\i>t ) ), (B.4) 
whose formal solutions are written as 

W ( L ) ) = \^)+Gf ) V\^)- (B.5) 

Although the Green's operators and the wave functions above are complex values, we do 
not necessarily have to handle complex values when the potential V{x) is real. For this we 
define the principal values of the two-body Green's operators VGl and VGo.l 

VGl = v e—t^TW) (R6) 

As is Gq 2, an analytical form of VG a 2 is known and these operators are related as 

G^l=VG ,L^iq^\jL}{j L \, (B.8) 
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A scattering wave function corresponding to VGo,l, namely standing wave solution \^>l) 
satisfies 

= \3l)+VG %l V\$ l ), (B.9) 
and a formal solution of this is given as 

\^l) = \3l)+VGlV\jl). (B.10) 
From the standing wave solution, the outgoing and the incoming solutions are obtained as 

\* ( t ) ) = i^crhM. (B.n) 

where ICl is the scattering jf-matrix defined by 



>Cl = -q^{jL\V\i> L ), (B.12) 



which becomes tan S with a phase shift parameter S. Using the relations above, one obtains a 
relation between G^ L ' and VGl as 



which reduces to Eq. (|B,8|I if V(x) was 0, leading to ipL{x) = jr<(qx) and KLl = 0. 

Next, we discuss about asymptotic form of the Green's functions. The asymptotic forms of 
G { ^1 and VG ,l are obtained from their analytical forms as 

G&^-q^\h¥}{jL\, (B.14) 
vn 

TG , L ^q^\n L )( jL \. (B.15) 



These equations and the resolvent equations together with the formal solutions Eqs. (|B,5|I and 
(|B~T0|) lead to 

Of 5 - ~q^\h^m\, (B.16) 

VG L ^q-^\n L ){i> L \. (B.17) 
ft 

Finally, we describe how to calculate Eq. (|2 . 23 p . which we write simply as 

V (x) = {x\G ( L +) \u;). (B.18) 



Using Eq. (TB7T3I . one can write rj(x) as 

Tj(a!) = fj(x) - ^^^(a:)-— L_<^|<&>, (B.19) 

ft 1 — lK,L 

where a new function fj(x) is defined by 

f)(x) = {x\VG L \u). (B.20) 
From Eq. (|B.17|I . the asymptotic form of f](x) can be written as 



7TI 

V{x) -* <7T2 n i( l ? a: )< t /'Lk) (B-21) 
x— >oo ft 



In actual calculation, the function r/(x) is obtained by solving the ordinary differential 
equation 

[E q -T L (x)-V(x)]fj(x)=u>(x) (B.22) 

with the boundary condition 

fj{x) oc nh{qx). (B.23) 

x — too 

These relations give the asymptotic form of rj{x) as 

^^oo^^I^W ( ^ } - (B - 24) 
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Appendix C: Mesh points for x and y variables 
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Crucial procedures in our numerical calculations are to solve the differential equations Eq. 
(|2.30[) and the Fourier-Bessel transformation Eq. (|2.28[) , which are related to x- and y-mesh 
points, respectively. In this appendix, we give some remarks on these mesh points. 

Both mesh points are taken in uneven distances so as to be shorter near the origin to take 
into account of short range nuclear potentials. 

Uneven mesh points, for x-mesh, e.g., are created with the same functional form as the one 
used in Ref. [5] 

c(x + to)x 



t(x) 



X + So 



or inversely 



B (t) = 



-(Ct ~t)- y/(ct ~ t) 2 + 4cS t 

2c ' 



(C.l) 



(C.2) 



with equidistant i-mesh points. The parameters of Eq. (|C.1|) . c, to, and so, are determined from 
the following conditions: 

1. The a;- mesh size near the origin : Axo 

2. The a;-mesh size at large distance (the infinity) : Ax x 

3. A value of a;-mesh points, x m , and the number of the mesh points for < x < x m : N x 
We can choose the distance of t-mesh points, At, as arbitrary, and thus 



N x At = t{x n 



From Eq. JOT 



dt 

dx 
dt 
dx 



cto 
so 



(C.3) 

(C.4) 
(C.5) 



Then, 



go 
ct 



At 



1 



= -At 

c 



(C.6) 
(C.7) 



Ax 

AXoo 

Using the values of N x , x m , Axo, and Axoo as an input, we rewrite the above conditions 

At 



Ax c 

N x - 



so 



to 



Mm fj 

A XQ JVa: 



AXqc 

Axo 



(C.8) 
(C.9) 
(CIO) 



In the present calculations, we set At = 1 (fm) both for x- and y-mesh points; N x — 60(100), 
x m = 10(80) fm, Axo = 0.025(0.033), Ax x = 0.3(1.25) for x- (y-) mesh points. 
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